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A grid of numerical simulations of double-diffusive convection is presented for the astrophysical case where viscosity 
(Prandtl number Pr) and solute diffusivity (Lewis number Le) are much smaller than the thermal diffusivity. As in 
laboratory and geophysical cases convection takes place in a layered form. The proper translation between subsonic flows 
in a stellar interior and an incompressible (Boussinesq) fluid is given, and the validity of the Boussinesq approximation 
for the semiconvection problem is checked by comparison with fully compressible simulations. The predictions of a 
simplified theory of mixing in semiconvection given in a companion paper are tested against the numerical results, 
and used to extrapolate these to astrophysical conditions. The predicted effective He-diffusion coefficient is nearly 
independent of the double-diffusive layering thickness d. For a fiducial main sequence model (15 Mq) the inferred 
mixing time scale is of the order 10 10 yr. An estimate for the secular increase of d during the semiconvective phase is 
given. It can potentially reach a significant fraction of a pressure scale height. 
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In models of stellar structure, situations are found where 
the heavier products of nuclear burning provide stability 
to a zone which otherwise would be unstable to convective 
overturning. Such a zone, or part of it, would become con- 
vective if something managed to mix its composition (R.J. 
Tayler 1953). The question whether such a zone should be 
treated as if it were mixed or not has become known as 
the semiconvection problem. Answers to this question dif- 
fer substantially. In practice, recipes are used containing a 
free parameter that allows the degree of mixing to be var- 
ied. Calculations in which such a parameter is adjusted to 
match observations are then called 'with semiconvection'. 
Commonly used prescriptions are those of Langer 1985 and 
Maeder 1997. 

The presence of a semiconvective zone has only a minor 
effect on the thermal structure of the star. The assumed 
amount of mixing of composition is critical, however, be- 
cause the evolution of the star is sensitive to the precise 
distribution of products of nuclear burning with depth in 
the star. The main goal of a theory for semiconvection is 
thus a good determination of the rate of mixing. From the 
perspective of the stellar evolutionist, the theory would ide- 
ally provide a formula for the rates of mixing and transport 
(the effective diffusivities), as functions of local thermody- 
namic state and composition, and their gradients. 

In Spruit (1992, hereafter S92) such formulas were 
derived, adapting the known physics of double-diffusive 
convection (Turner 1979, 1985, Proctor 1981, Huppert & 
Turner 1981, Schmitt 1994) to the case of a stellar interior 
(discussed first in this context by Spiegel 1969, 1972). The 
expression developed in S92 makes use of simplifications 
valid in the limiting case of very large Rayleigh number and 
very low solute diffusivity. In a companion paper (hereafter 



SI 3), the analysis in S92 is extended to cover the more 
moderate conditions accessible with numerical simulations. 
In the following, mixing rates and their dependence on the 
intrinsic parameters governing the problem are measured 
with such simulations. The predictions of S13 are tested 
against these results, and then used to estimate the ex- 
pected mixing rates in semiconvective zones of stars. 

One of the predictions in S13 is the existence of a max- 
imum density ratio (the ratio R p of stabilizing solute to 
destabilizing thermal stratification) for which a steady lay- 
ered state is possible. In a slightly different guise, this limit 
also figures prominently in Proctor's (1981) analytic analy- 
sis. In this analysis, he proved that in the limit of vanishing 
solute diffusivity there exists a layered state at any Rayleigh 
number above the critical value for convection in the ab- 
sence of a stabilizing solute gradient, provided the density 
ratio is below this critical value. The model in S13 does not 
predict what happens for density ratios just above or be- 
low this maximum value. The numerical results presented 



in Sect. 5.2 clarify how the system behaves in this case. 



The development of a linear initial gradient into the fi- 
nal state of overturning layers separated by diffusive steps 
is studied with a few examples. This transient process 
shows the 'Kato-oscillations' expected from linear theory 
(see Sect. 5 for an example). It is not very relevant for as- 
trophysical application, however, since the transition to the 
layered state happens on a time scale (a finite number of 
buoyancy oscillations) that is negligible compared with the 
time scales of interest, and is bound to depend on the de- 
tails of the initial state. Instead, the focus here is on the 
transport of heat and solute in a double diffusive staircase 
with layers of given thickness, as a function of the intrinsic 
parameters of the problem. 
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In Sect. 2 the known physics of double diffusive con- 
vection of the semiconvective type is reviewed, in general 
and from an astrophysical perspective. Sect. 3 describes 
the transport properties of a single double-diffusive layer in 
terms of the model of S13. The numerical methods are given 
in Sect. 4. In Sect. 5 the results of our parameter study are 
shown and compared with the predictions of this model. 
Application to the case of semiconvection in a 15 Mq main 
sequence star is discussed in Sect. 6. 

2. Semiconvection and double-diffusive convection 

Situations where a fluid is stabilized by the density gradi- 
ent due to a dissolved heavy constituent occur in nature. 
An example is convection in arctic oceans (fresh meltwa- 
ter cooled from above, stabilized by the salts dissolved in 
the sea water, see the review in Schmitt 1994). Intensively 
studied are East-African rift lakes (lakes Kivu, Nyos and 
Mormon, cf. Schmid et al. 2010). These are heated from 
below by volcanic activity, which also is a source of dis- 
solved gases (carbon dioxide and methane, hereafter the 
'solute'). Their density stratification is stabilized against 
convection by the stable gradient resulting from the weight 
of the carbon dioxide. Efforts to prevent catastrophic re- 
lease of carbon dioxide (Lake Nyos, e.g. Sigvaldason 1989) 
or to enable safe commercial exploitation of methane (Lake 
Kivu, Nayar 2009) have led to extensive study of the fluid 
flows, heat flux and mixing rates in these natural double- 
diffusive systems. 

The gradients in temperature and solute in the East- 
African lakes and the arctic are observed to be 'stepped': 
consisting of a stack (called 'staircase') of thin layers 
(decimeters to decameters). Inside a layer, overturning con- 
vection keeps the composition nearly uniform, with sta- 
ble gradients in temperature and composition separating 
the layers. The physics involved is easily reproduced under 
controlled laboratory conditions (Turner & Stommel 1964, 
Turner 1985) Q The layers are very long-lived: of the order 
of months or more in the geophysical examples mentioned, 
orders of magnitude longer than the convective turnover 
times inside the layers. 

In the stable gradients between the overturning layers 
the transport of the stabilizing solute takes place by diffu- 
sion instead of convection. This strongfly limits the effective 
transport of solute through the double-diffusive staircase. 
Residence times on the order of 1000 yrs are inferred for 
the solutes in lake Kivu, for example (Schmid et al. 2010). 
This is 8 orders of magnitude longer than the convective 
turnover times in these layers. The transport of heat is also 
strongly reduced; this is exploited for heat storage in solar 
ponds (cf. Lu and Swift, 2001). Similarly low fluxes have 
been measured in thermohaline staircases in the arctic and 
antarctic oceans (e.g. Padman & Dillon 1987). 

Theoretically, the observed layered nature of double- 
diffusive systems is well understood. Central to this un- 
derstanding is the fact that linear stability analysis does 



1 Also on a coffee table. Latte macchiato in a tall glass often 
shows the effect nicely. After the coffee is added to the milk, a 
stably stratified gradient of milk/coffee mix develops (showing 
internal gravity waves in the form of a sloshing motion with 
a period of a few seconds). After about a minute, the initially 
smooth gradient starts dividing into thin (a few mm) layers, 
visible at low contrast. In the course of several minutes these 
merge into a smaller number of more clearly defined layers. 



not provide relevant clues to their behavior, because the 
double-diffusive case of thermal convection stabilized by a 
slowly-diffusing solute is subcritical. That is, a nonlinear 
form of the instability, in the form of a stable overturning 
flow, exists already below the critical temperature gradient 
for onset of linear instability. The linear stability condi- 
tion is thus not relevant for the behavior of the system (cf. 
Schladow et al. 1987). 

Linear instability predicts internal gravity modes to set 
in above some critical value of the temperature gradient, 
growing in amplitude by the effect of thermal diffusion: the 
so-called Kato oscillations (Kato 1966). Such oscillations 
(cf . movie at Fig. [5 ) transport a negligible amount of heat 
or solute, compared with overturning motions of the same 
amplitude. For this reason alone, linear stability arguments 
cannot be used for estimates of the mixing rate in semi- 
convective zones. More important, however, is the subcrit- 
ical nature of double-diffusive convection. Proctor (1981) 
shows analytically that, in the limit of vanishing diffusivity 
of the solute, the layered form of convection exists whenever 
the Rayleigh number exceeds the critical value for ordinary 
convection, irrespective of the strength of the stabilizing 
component. This assumption of vanishing solute diffusivity 
is eminently satisfied in the astrophysical case and holds 
reasonably well in geophysical and laboratory experiments. 

In Stevenson 1979 it was assumed that the nonlinear de- 
velopment of the overstable oscillations would lead to satu- 
ration of the wave at a finite amplitude. This has been the 
rationale for some prescriptions used in astrophysics (e.g. 
Langer et al. 1985). The assumption of saturation at finite 
amplitude is appropriate for the more common supercritical 
forms of oscillatory instability, but not for the subcritical 
case, where a finite amplitude state exists for parameters 
where the system is still linearly stable. The assumption 
was also in conflict with observations: the laboratory and 
geophysical systems all showed the same characteristic form 
of convection in a system of overturning layers (e.g. Turner 
& Stommel 1964), none that settled into finite amplitude 
oscillations. 

In the literature on semiconvection it is often argued 
that the Prandtl number in astrophysics is much lower than 
in the geophysical and laboratory cases. The implied as- 
sumption that the understanding of double diffusive con- 
vection developed in these contexts can be set aside, is 
not necessary however. Proctor's analysis, for example, is 
largely independent of Prandtl number. It is valid as long 
as Pr is not larger than of order unity, and viscosity is 
larger than the solute diffusivity. This is also satisfied in 
the astrophysical case. 

In physical terms, the reason for the subcritical behavior 
can be understood with an energy consideration (S92). The 
amount of energy it takes to overturn a layer of thickness d 
against a stable gradient scales as d 2 (as in a harmonic oscil- 
lator at the buoyancy frequency of the stratification). The 
expense in initial energy needed per unit mass to put the 
system into its finite-amplitude, layered state thus vanishes 
as d, down to some value where damping losses stabilize 
the system. A small initial perturbation, or an initial Kato 
oscillation, is sufficient to provide the energy for overturn- 
ing into thin layers. Once established, this layered state is 
a stable form of convection. In fact, reproduction of the os- 
cillatory phase in the laboratory requires very careful setup 
of the initial stratification (Shirtcliffe 1967, see also the dis- 
cussion in Huppert & Turner 1981). 
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This agrees with the observation in laboratory exper- 
iments and geophysical systems like lake Kivu mentioned 
above, i.e. that the layering first sets in at a small thickness 
(cf. footnote above). In this context, the formation of layers 
found in the simulations by Rosenblum et al. 2011 is in line 
with expectation. Much less well-defined is the evolution on 
longer time scales, in particular the question how and on 
which time scales the layer thickness evolves (see Sect. 6.6 
and S13). 



2.1. Double diffusive convection in stars 

Double diffusive convection in stars has traditionally been 
regarded as a piece of physics to be treated separately 
from the geophysical examples, since the numerical val- 
ues of controlling parameters such as the Prandtl number 
are quite different. Apart from the difference in equation 
of state, however, the hydrodynamic equations are identi- 
cal. Differences in physics that might be present between 
the two cases are in fact not apparent in the elementary 
recipes for semiconvection used in stellar evolution codes. 
The failure of these recipes when applied to the geophysi- 
cal case is traditionally not considered an argument against 
their application in stars. Such recipes include (a) to assume 
that no mixing takes place at all (the 'Ledoux' recipe), an- 
other (b) to ignore the stabilizing effect of the solute gra- 
dient ('Schwarzschild' recipe, yielding a very high mixing 
rate) (c) some interpolation between these recipes, (d) to 
assume that the amount of mixing is such that the layer be- 
comes marginally stable to overturning, and (e) the some- 
what more physically motivated oscillation-based recipe of 
L anger et al. above. 

Define the Prandtl number Pr = v / 'kt and the Lewis 
number Le = ks/kt, where kt is the thermal diffusivity 
and Kg the diffusivity of the stabilizing 'solute' (He diffusing 
in H, say). Because of the high thermal diffusivity mediated 
by photons, Pr and Le are very small numbers, some 8-10 
orders of magnitude less than in geophysical cases. 

Such small parameter values cannot be covered realis- 
tically in numerical simulations. If t c is a typical convec- 
tive time scale (as estimated from the superadiabaticity 
and pressure scale height in the semiconvective zone), very 
small length scales, of the order (t^ks) 1 ^ 2 would have to 
be resolved to represent the interaction between flow and 
diffusion. This is not computationally possible at present 
even in 2-dimcnsional simulations. 

Translation from the numerically accessible parameter 
values to an astrophysically relevant parameter range there- 
fore requires scaling of the results over the orders of magni- 
tude in-between. A valid extrapolation cannot be found by 
mere intuitive inspection of numerical results, since these 
are too far from the target regime. On their own, the ap- 
plicability of numerical results to the astrophysical case is 
bound to remain diffuse, or applicability as an explicit goal 
has to be given up altogether. 

Some physical understanding that includes this asymp- 
totic regime is necessary. A theory that accomplishes this 
is given in the companion paper (S13). This is made possi- 
ble by making explicit use of the observed separation into 
layers of convective overturning between stable diffusive in- 
terfaces. 

An important difference between the astrophysical and 
geophysical cases is the equation of state. For the thin lay- 



ering expected, sound travel times are very short compared 
with convective time scales. As a consequence, the fluid in 
a semiconvective zone behaves as nearly incompressible. A 
Boussinesq approximation can then be used for the calcula- 
tions, provided a small complication is properly taken into 
account. Whereas in the incompressible case the convective 
and diffusive heat fluxes are both governed by the tem- 
perature gradient, convection and radiative diffusion are 
governed by different quantities in the compressible case 
(entropy gradient and temperature gradient). This affects 
in particular the definition of the Nusselt number as a mea- 
sure of the efficiency of heat transport. An exact translation 
between the se ca ses is possible (Massaguer and Zahn 1980, 
see below in 6.1 ). 

A second difference concerns the status of the heat flux 
in the formulation of the problem. The effect of semiconvec- 
tive mixing on the star's structure in evolution calculations 
is found to be small during the semiconvective phase itself, 
somewhat independent of the way semiconvection is ap- 
proximated. (This is in part because of the limited extent 
of a semiconvective zone). Its effect on the radial profile 
of elemental composition, however, is of lasting importance 
for later evolutionary stages (see Langer et al. 1989). 

The consequence is that, in contrast with laboratory 
and geophysical situations, in a stellar model the heat flux 
F, rather than the temperature gradient can be consid- 
ered as known. Since the radiative contribution F r to the 
heat flux is known to good approximation from the thermal 
structure of the star, the convective heat flux F c — F — F x 
transported by semiconvection is also known. The efficiency 
of convection: i.e. how close the mean thermal gradient is 
to the adiabatic gradient, follows from the imposed heat 
flux (instead of the other way around as in a laboratory 
experiment). We return to this distinction in Sect. 6.4. 

2.2. Size of the parameter space 

The fluid is described by the thermodynamic variables 
defining its local state (e.g. temperature and density) and 
the material functions (e.g. pressure, diffusivities, viscos- 
ity) . In addition the gradients of the thermodynamic vari- 
ables with depth, and the acceleration of gravity are rele- 
vant for the properties of the flow. Taken together, these 
quantities form a large parameter space, and it might be 
concluded that realistic numerical simulations of semicon- 
vection would have to be done for individual zones in indi- 
vidual stellar models. 

The equations of fluid dynamics have symmetries, how- 
ever, so that the independent degrees of freedom are far 
fewer. They can be represented by five dimensionless pa- 
rameters: a Rayleigh number Ra, the layer thickness d in 
units of the pressure scale height H , the Prandtl number, 
the Lewis number, and a density ratio R p which measures 
the ratio of the stabilizing (solute) gradient to the desta- 
bilizing thermal gradient. The behavior of semiconvection 
at some point in a star can be defined in terms of these 
parameters. 

The Boussinesq approximation corresponds to the limit 
e = d/H <C 1. The pressure scale height disappears as 
parameter in this limit; all dependence on d is subsumed in 
Ra. This reduces the number of parameters of the problem 
to four. 

By a fortunate coincidence, it turns out that as long as 
the Prandtl number is less than unity, the results are effec- 
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tively independent of Pr. This further reduces the number 
of independent parameters to only three. Since measure- 
ment of the mixing rate in each individual case does not re- 
quire a very expensive simulation, this allows a significant 
volume of parameter space to be covered, and a comparison 
with the model predictions in S13 to be made. 



3. The layered state 

3.1. Layer formation, layer thickness 

The formation of a layer from an initially smooth and static 
gradient starts with a well-known oscillatory instability, 
with oscillation periods of order the buoyancy period of the 
stratification. This is followed by nonlinear development 
into an overturning flow. The transition can be observed in 
numerical simu lations (Merryfield 1995, Rosenblum et al. 
2011, Sect. |5.5| below), and in very carefully designed labo- 
ratory experiments, but is not seen in geophysical cases. An 
unperturbed smooth initial gradient in solute and temper- 
ature is an artificial case that is unlikely to be realized in 
nature. In addition, the intermediate oscillatory state lasts 
for only a finite number of oscillation periods (5-10 in the 
results below), a very short time scale compared with the 
life time of the double diffusive layers. In the case of lake 
Kivu, for example, buoyancy periods are in the range 5-30 
min, the life time of individual layers a few months. 

In laboratory and kitchen table experiments the lay- 
ers are initially very thin. This is understood from the en- 
ergy argument above. The layers slowly merge, either by 
fading of contrast between neighboring layers, or the ver- 
tical migration of interfaces towards adjacent layers. The 
details of this process have not been studied much (but 
see McDougall 1981, Young & Rosner 2000, Ross & Lavery 
2009). A plausible estimate of the rate of growth of the 
layer thickness can be given in terms of the effecti ve s olute 
diffusivity of the system, however (S13, and Sect. 6.5). 



In the simulations reported here, the evolution of layer 
thickness is not included. The thickness (d) is therefore 
treated as a free parameter of the problem. It enters 
thro ugh the Raylcigh number (in the Boussincsq limit, see 
2.2 above). In the astrophysical application of imposed heat 
flux, however, the resulting mixing rate is effe ctive ly inde- 
pendent of the layer thickness (S13 and Sect. 6.5 below). 
This is a fortunate circumstance, since following the evolu- 
tion of a stack of layers through its merging processes would 
require far lengthier simulations (for an example see Young 
& Rosner 2000). 

3.2. Structure of a layer 

An individual layer in the double diffusive staircase con- 
sists of a zone of overturning convection, separated from 
its neighbors by stably stratified stagnant zones. In a stag- 
nant zone the transport both of heat and solute takes place 
by diffusion, hence the profiles of S and T are steep and 
approximately linear with depth in this zone. In the over- 
turning zone the gradients are shallow, except in the thin 
boundary layers at its interface with the stagnant zone. The 
thickness d s of the stagnant zone is usually only a small frac- 
tion of the layer thickness d, but can be as high as 20% close 
to marginal conditions for layer formation (for a discussion 
of its dependence on parameters see S13). 



Apart from deformations by internal gravity waves in 
the stagnant zone, both the vertical and the horizontal ve- 
locity components vanish at its interface with the overturn- 
ing zone. Apart from the small amount of solute carried 
by the flow, the overturning zone thus behaves essentially 
like laboratory convection in a box, the interfaces with the 
stagnant zones acting like the top and bottom plates, with 
viscosity enforcing no-slip conditions. Because of the peri- 
odic boundary conditions in the horizontal direction, free 
slip conditions would allow the stagnant zones as a whole to 
start moving sideways. To the overturning flow, the stag- 
nant zones would still present internal no-slip interfaces, 
however. We have chosen the no-slip conditions since the 
freedom of such large-scale flows is somewhat of an arte- 
fact of the periodic boundary conditions in the horizontal 
direction. In layers of realistically large horizontal extent 
they would probably not happen. This picture has already 
been put forward early on in the interpretation of labora- 
tory and geophysical observations (Shirtcliffc 1967, Linden 
& Shirtcliffe 1978) and of the results from numerical simula- 
tions such as those reported below. It has become the stan- 
dard interpretation used for deducing heat and solute fluxes 
in thermohaline staircases in geophysics (e.g. Padman & 
Dillon 1987, Turner & Stommel 1964, Schmid et al. 2010). 



3.3. The overturning zone 

At the boundary with the stagnant zone the overturning 
flow has three nested boundary layers: for temperature, so- 
lute, and flow speed; the thermal, solute and viscous bound- 
ary layers respectively. Their thicknesses are determined by 
the thermal diffusivity kt, solute diffusivity ks, and (kine- 
matic) viscosity v. The highest diffusivity (thermal) has 
the widest boundary layer. For astrophysical conditions, 
kt 3> Kg, v 1 so the solute and viscous sublayers are thin 
compared with the thermal boundary layer. 

The horizontal velocity vanishes in the stagnant zone 
(apart from internal wave modes), so it imposes no-slip con- 
ditions on the overturning flow. Since the viscous bound- 
ary layer is thin compared with the thermal boundary 
layer, however, this boundary condition has little effect on 
turnover times and the resulting heat flux; it becomes no- 
ticeable only at low Raylcigh numbers. 

The need to properly resolve the narrowest of the 
boundary layers in a numerical simulation determines the 
number of grid points needed in the vertical direction. In 
the astrophysical case, viscosity is typically of the same 
order or somewhat larger than solute diffusivity, so the re- 
quire d re solution is set by the solute boundary layer (cf. 
Sect.Ol. 



3.4. Heat transport 

The flow in the overturning zone is driven entirely by 
boundary layers at the top and bottom steps of the layer, 
much like in laboratory convection in a box (Niemela et 
al. 2000). Except for these thin boundary layers, the (hor- 
izontal averages of) entropy and composition are almost 
uniform inside the layers. Under the conditions in a stel- 
lar interior, the thermal diffusivity is much larger than the 
diffusivity of the solute (Helium in Hydrogen, say). The 
thickness of the solute boundary layers is then much smaller 
than the thermal boundary layers, and the amount of so- 
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lute in transit across the layer vanishes in this limit. The 
convective flux is then almost the same as in the absence of 
a stabilizing solute (cf. Schmitt 1994). Under these asymp- 
totic conditions, the dependence of heat flux on Rayleigh 
number can be taken from a simple estimate, as was done 
in S92, or a result from laboratory measurements of con- 
vection can be used (e.g. Niemela et al. 2000). With the 
convective flow thus known, the flux of solute can be cal- 
culated as well. 

Convection experiments in gaseous Helium at tempera- 
tures just above the critical point by Castaing et al. 1989 
yielded the fit 

Nu T = 0.23 Ra 282 . (1) 

over the range 10 7 < Ra < 10 14 . More recent measure- 
ments by Niemela et al. 2000 using superfluid 3 He, over 
the remarkable range 10 6 < Ra < 10 17 gave a marginally 
different result: 



Nu T = 0.124Ra' 



0.309 



(2) 



This can be compared with the estimate based on a 2- 
dimensional argument in S92: 



Nu T » 0.5 Ra' 



0.25 



(3) 



Since the Prandtl number in the laboratory experiments (of 
order 0.7) is not far from unity, Ra in (jTJ [2]) can be identified 
with Ra* in Expression |3| then agrees within a factor 
2 with (||), up to Ra* = 10 15 . 

The expressions above are of the form Nut = aRa b . 
They are fits for large Ra; to better cover lower Rayleigh 
numbers as well, we modify the expression slightly: 



Nu T - 1 = a (Ra* - Ra* c ) 6 



(4) 



where Ra* c = PrRa c , and Ra c ~ 1400 the critical Rayleigh 
number for convection with no-slip boundary conditions. 
With this change the Nusselt number approaches the dif- 
fusive value of 1 as Ra approaches the critical value. 

3.5, Solute transport 

The transport of solute is governed by the diffusion from 
the overturning flow into the steep gradient in the stag- 
nant zone. This happens in a narrow boundary layer of the 
flow. Its thickness 5 S is of the order <5g = (kst) 1 / 2 , where 
r is the time the flow travels in contact with the bound- 
ary before descending/ascending back into the interior. In 
the same way, the thermal boundary layer has thickness 
St = (ktt) 1 ' 2 - The fluxes of solute and heat carried in 
these boundary layers also determines the flux across the 
overturning zone as a whole. This predicts that the solute 
and heat fluxes Fs, Ft are related by 



F S /F T ~ (k s /k t ) 1/2 . 



This is the well-known relation derived in various ways for 
double-diffusive convection (cf. Turner 1985). In a more ac- 
curate form it appears in the model of S13 and in the nu- 
merical experiments reported below. 

The amount of solute transported across the overturn- 
ing zone is limited by the fact that convective plumes de- 
velop only from material with a net buoyancy of the unsta- 
ble sign. This limits the density contrast of the stable solute 
that can be carried to a fraction 1/ R p of the density con- 
trast across the layer as a whole (Schmitt 1994, S92, S13). 



Solute with a higher density contrast is not mixed into the 
convective flow: it remains in the stagnant zone. The solute 
flux carried by convection in the overturning zone has to 
match the diffusive flux in the stagnant zone, and idem for 
the heat flux. In S13 it was shown that these conditions, 
together with expression Q form a complete model, de- 
termining the effective transport coefficients as well as the 
thickness of the stagnant zone, as functions of the three 
parameters of the Boussinesq problem, R p , Ra*, Le. 

In the estimates above the amount of solute transported 
is given by the amount flowing through the solute boundary 
layer of thickness 5 S . This assumes that the interaction with 
the stagnant zone is determined entirely by diffusion of the 
solute. This is a minimum: the stagnant zone as a whole is 
stable in the sense of the Richardson condition (one verifies 
that this is equivalent to the fact that R p > 1), but near its 
boundary with the overturning zone the density gradient is 
lower, and a certain degree of mixing by shear instabilities 
is to be expected. This will increase the amount of solute 
that is of sufficiently low buoyancy to be carried with the 
flow. As argued in S13, this additional amount scales with 
5 S , so that its net effect is equivalent to an increase of 5 S 
by a factor q, of order unity, increasing the solute flux by 
the same factor. The predicted relation between solute and 
thermal Nusselt numbers derived in S13, with this erosion 
factor included is 



Nus - 1 = 



Le 1 / 2 ;?, 



(Nu T - 1) (R P < Le 



-1/2N 



(6) 



The value of q will be a fitting parameter in the quantitative 
comparison with the numerical results. [The value R p = 

1/9 

Le ' is the largest for which the model predicts existence 
of a layered state at any Rayleigh number, see also Sect. 
5.2 below. 1 



3.6. Model summary 

The model used for comparison with the numerical simula- 
tions is defined uniquely by the ingredients described above. 
It makes use of only minimal assumptions about of a double 
diffusive layer: (a) convection in the overturning zone of a 
layer is described by a fit to laboratory measurements, (b) 
transport in the stagnant zone is by diffusion only, and (c) 
the fluxes of solute and heat are related by the buoyancy 
limit in the overturning zone (S13 eq. 17). These together 
determine the thickness of the stagnant zone and the effec- 
tive transport coefficients of heat and solute through the 
layer as functions of the parameters of the problem, Le, Ra 
and R p . The behavior of the model is discussed in more 
detail in S13. 



(5) 4. Numerical simulations 



All equations were calculated on a 2D rectilinear Cartesian 
grid in terms of finite differences. The horizontal coordi- 
nate is labeled with x, the vertical coordinate with z, where 
z — is the bottom and z = 1 the top boundary. The set 
of equations are implemented into the ANTARES software 
framework (Muthsam et al. 2010). Advective currents are 
treated by a weighted essentially non-oscillatory finite vol- 
ume scheme in fifth order (Shu & Osher 1988), the phys- 
ical diffusion is handled by a fourth-order finite difference 
discretization. A second order total variation diminishing 
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scheme is chosen as time integrator. To avoid odd-even de- 
coupling, a MAC grid, which locates vector variables at cell 
faces and scalar variables at cell centres, is used. 

For a detailed description of the numerical solution of 
the binary mixture equations presented here see Zaussinger 
2011. 



4.1. Boussinesq approximation for a binary mixture 

In the Boussinesq approximation the continuity equation is 
reduced to that of an incompressible fluid, 



dp 
dt 



= 0, 



(7) 



such that in the equation of motion the density is taken to 
be a constant po, 



du 

~dt 



1 

Po 



VP + (-a T Q + a s S)g + V ■ {vVu), (8) 



where the 'expansion coefficients' ax, as describe the den- 
sity effects of variations O, S in temperature and solute, as- 
sumed small compared with the mean temperature and 
salinity S. They are given by advection-diffusion equations: 



d9 

dT 

dS 
di 



= —u ■ VO + V • (k t V9), 



= -m-V5 + V-(k s VS*), 



(9) 



(10) 



where v is the kinematic viscosity. The mean gradients V0, 
V/S can be expressed more usefully in terms of the buoyancy 
frequencies: 

« T g'Ve, (ii) 



A s 2 = a s g • VS, 



(12) 



The density ratio R p is then defined as in the compressible 
case, eq. (32). 



The numerical algorithm solving this set of equations is 
based on a semi-implicit scheme. Intermediate values for the 
velocity field u* are calculated explicitly from the equations 
of motion. By the nature of the incompressible equations, 
the pressure update is done implicitly, by solving a Poisson 
equation: 



AP 



At 



(v-o 



(13) 



The resulting pressure P leads to the required diver- 
gence free velocity field at the new time step n + 1. 



At. 

Po 



VP. 



(14) 



4.2. Compressible fluid equations for a binary mixture 

Verification of the Boussinesq results has been done with 
simulations of the fully explicit compressible fluid equa- 
tions. The fluid is assumed to be an ideal gas, which is a 
good approximation to a binary gas mixture of our interest. 
These are the continuity equation 



dp 
di 



+ V • (pu) = 



(15) 



the partial density equation 
d(pc) 



dt 



V • (pen) = V • (pK c Vc) 



(16) 



the momentum equation 
d{pu) 



dt 



V • (puu + FT) = pg + V • r 



(17) 



the total energy equation 
dpE 



dt 



V-[u(pE + P)} =V-(K h VT) + V-{uT) + pgu (18) 



and the equation of state 



P : 



KpT 

p(c) 



(19) 



where c is the solute mass fraction, E — \\u\ 2 + e is the 
total specific energy in units of energy per mass, e is the 
specific internal energy, is the heat conduction coeffi- 
cient, k c is the solute diffusion coefficient, r is the viscous 
stress tensor, g the gravitational acceleration, I the unit 
tensor, fi = p(c) the mean molecular weight, and 1Z the gas 
constant. 



4.3. Units, boundary and initial conditions 

As unit of length we use, for the Boussinesq cases the layer 
thickness d, for the compressible calculations the pressure 
scale height H. A nominal convective turnover time is used 
as unit of time. As units of temperature (potential tem- 
perature in the compressible case) we use 1/cxt, for so- 
lute concentration I/as- The density ratio then becomes 
R p = V5/VT. The Boussinesq equations are invariant to 
arbitrary additive constants in temperature and solute. We 
set these such that T and S are zero at the top boundary. 

Because of the symmetries of the problem, there are 
fewer independent parameters than physical variables de- 
scribing it. Hence some of the physical quantities appearing 
in the problem can be set to unity. We choose for these: 
the temperature difference between top and bottom of the 
layer, the density at the bottom of the layer and the accel- 
eration of gravity. The Rayleigh number Ra* is then con- 
trolled through the thermal diffusivity kt, the solute differ- 
ence between top and bottom through the density ratio R p , 
and the solute diffusivity through the Lewis number Le. 

Most of the calculations were done in a box simulat- 
ing a single layer from the double-diffusive staircase, so the 
top and bottom boundaries coincide with the steps between 
layers. This ignores the distortions of the interfaces by sur- 
face waves, but since the essence of the double layering 
phenomenon is that the transport across the interface is by 
diffusion, this is not expected to make a big difference. To 
check that this is indeed the smaller set of simula- 

tions was done in which a step is present inside the volume 
(Sect.[5~7|. The vertical boundary conditions are thus taken 
to be impermeable and stress-free. In the horizontal direc- 
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tion periodic conditions are used. 





= 


(z = 


0,1) 


(20) 


du x 






0,1) 


(21) 


dz 


= 


(z = 


S 


= Rp 


(z = 


0) 


(22) 


T 


= 1 


(z = 


0) 


(23) 


S 


= 


(z = 


1) 


(24) 


T 


= 


(z = 


1) 


(25) 



As initial condition the stratification of temperature and 
solute was taken to be horizontally uniform with cither 
a linear gradient between the values at top and bottom 
(the 'linear' case below), or something approximating the 
boundary layer structure expected of the final state (the 
'step' case). Small initial random perturbations are applied 
on the solute field. 

The numerical algorithm for setting up the initial con- 
ditions for the compressible fluid equations is an extension 
of the procedure presented in Muthsam et al. (1995, 1999). 

4.4, Numerical setup 

Being the thinnest boundary layer in the problem, 5g deter- 
mines the numerical resolution needed near the boundaries. 
Since the fine structure in the interior of the layer largely 
consist of boundary layers 'peeled off' from the boundaries, 
the same resolution is needed in the interior of the layer 
as well, and there is no need or justification for using non- 
uniform grids. 

The expected solute boundary thickness for Ra* = 
1.6 x 10 5 and a Lewis number of Le = 0.1 is 1.6% of the 
layer thickness d. With the high-order spatial discretiza- 
tion used, a nominal resolution of tt-bl = 3 points is needed 
across the critical structures that must be captured. The 
transport of solute across the layer is determined by the 
vertical structure of the boundary layers. Resolving these 
with 7Jbl = 3 translates to n z = 200 points needed in the 
vertical direction for this case. 

To test the actual convergence of the results with re- 
spect to resolution, a series of simulations with n z — 100, 
200, 300 and 600 has been run for Ra* = 1.6x 10 5 , Le = 0.1, 
R p = 1.15, see Table m Conver gence to 15% is reached al- 
ready at n z — 100. For the lowest Lewis number used in 
the results reported below in Table 2, Le=0.01, the num- 
ber of grid points needed would be VT0 times higher. The 
number of points across the height of the box used in these 
simulations, n z = 300, is deemed sufficient for the level of 
accuracy we are aiming for. 



n 2 


Nu s 


Nu T 


100 


26 


9.5 


200 


27 


10 


300 


30 


10 


600 


30 


10 



Table 1. A test series shows convergence of the Nusselt 
numbers with respect to resolution (number of grid points 
n z ) across the height of the box. 



Most of the calculations were done with a horizontal- 
to- vertical aspect ratio of 2:1. A few tests with different 




Fig. 1. Flow structure in a double diffusive layer. The 
temperature field (top) is more diffuse than the solute 
('Helium') concentration, as a result of the high thermal 
diffusivity. 



ratios showed that this choice does not affect the measured 



Nusselt numbers significantly (Sect. 5.4) 



5. Results of numerical simulations 

Over a wide range in parameter values about 100 numerical 
simulations have been done in the Boussinesq approxima- 
tion and about 20 simulations have been performed with 
the fully compressible code. Only the cases closest to the 
parameter range of astrophysical interest (about 60) are 
reported here (Figs. [2]|4j Table 1). 

An example of the flow structure is shown in Fig.JT] with 
Pr = 0.1, Ra, = 5 x 10 5 , Le = 0.01, R p =1.15. It shows 
the key characteristic of double diffusive convection: the 
boundary layers and the plumes of the solute are narrower 
than those of temperature, on account of the low solute 
diffusivity. 



5.1. Dependence on Ra* and R p 

Fig. [2] shows the dependence of the measured Nusselt num- 
bers Nut and Nus on the parameters Ra* and R p , for the 
case Le= 0.01, Pr=0.1. The results show some fluctuation, 
as expected from the limited number of overturning times 
for which the simulations were run. This results in uncer- 
tainties in the Nusselt numbers of the order of a factor 
1.5. For comparison with the theoretical model predictions 
(dashed lines) the adjustable erosion factor q in eq. Q has 
been set to 1.5. Within the scatter, the predictions appear 
consistent with the numerical results to a factor of 2 or less. 
In fact, leaving out the erosion factor (i.e. setting q = 1) 
does not make the fit a lot worse. 

The upper right panel shows Nu s vs. Nut, where each 
curve represents the results for a fixed density ratio R p 
and concurrently increasing Ra* . The approximately linear 
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relation between the two shows that the ratio of thermal to 
solute transport depends primarily on Lewis number and 
density ratio, but not much on Ra, as predicted by The 
actual transport efficiency of both depends of course on the 
Rayleigh number Ra* as well (upper left panel). 

5.2. Behavior near the maximum density ratio 

The data in Fig. [2] show how the Nusselt number Nut de- 
clines with increasing density ratio. The theoretical analysis 
in S13 predicts a maximum density ratio R pmlix for the ex- 
istence of a (statistically) stationary layered state. It is a 
slowly increasing function of Ra* , asymptotically approach- 
ing the value Le ' . To test this a series of simulations was 
done at Ra* = 10 6 , Le= 0.1. The predicted value of i? pmax 
for this combination is 1.2 (the asymptotic value would be 
1.3). Fig. [3] shows the resulting Nusselt numbers as func- 
tions of time. 

For R p = 1.2 and 1.3, the heat flux in the simulations 
quickly settles to a statistically steady value after a tran- 
sient due to the initial conditions used. For the three highest 
values of R p , the initial state smoothly settles to the purely 
diffusive state Nt = 1, the faster and smoother, the higher 
the density ratio. Values in-between (1.5 and 1.8) are char- 
acterized by large fluctuations superposed on a slow general 
decline, as if the system is undecided between settling on an 
overturning or a static diffusive state. A precise boundary 
for the existence of an overturning state is therefore some- 
what hard to define: it appears to depend on the length of 
time over which the flow is followed. 



The time dependences of the Nusselt number shown 
in Fig. [3] suggest a value i? pma x ~ 1-4, somewhat larger 
than the theoretical value of 1.2. For density ratios slightly 
larger than 1.2, it looks at first as if a statistically steady 
state has been reached, but on a longer time scale a con- 
tinuing decline of Nut is observed. The behavior with in- 
creasing R p is therefore somewhat gradual, at least for the 
finite length of time over which the simulations have fol- 
lowed the development of the flow. Interestingly, the dis- 
tinction between diffusive and overturning final states be- 
comes sharper with increasing length of time over which the 
simulation is followed. This also affects the results in Fig. [2} 
At the largest density ratios shown there, the numerically 
measured Nusselt numbers are systematically higher than 
predicted. This can be attributed to the fact that at these 
R p the simulations were not run as long as those of Fig. [3] 

The numerical results thus confirm the predicted exis- 
tence of a maximum density ratio for the overturning lay- 
ered state, but with the added twist of very long settling 
times when R p is near the maximum value. This also an- 
swers a question that was left undecided in S13: the the- 
ory does not predict what happens close to the maximum 
density ratio. The results in Fig. [3] indicate that the sys- 
tem vacillates between overturning and diffusive states for 
increasingly long periods of time as the maximum is ap- 
proached. 

The fast settlement towards a diffusive state for high 
values of R p is not surprising given that the linear sta- 
bility criterion predicts exactly such a result for R p > 
(1 + Pr)/(Le + Pr) (e.g Veronis 1968, Stevenson 1979). 
For the parameter values of Fig. 3 this yields R p <; 1.8. 
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100 200 300 400 500 600 
N T . t 



Fig. 3. Time dependence of the thermal Nusselt number 
Nut in simulations with Pr=l, Le=0.1, Ra* = 10 6 , for 
density ratios (labeled) near the maximum value for which 
a layered state is predicted by the theoretical model. Time 
t in units of the thermal buoyancy time scale N^ 1 . 



The difference between 1.8, where Nut ~ 1 drops only ini- 
tially, but then only mildly so if at all, and a value of 2.0, 
which shows a continuous decay seems consistent with this 
classical result. 

The results also provide an interesting link to Proctor's 
(1981) analysis. In this theory the maximum value R p = 

1/9 

Le ' appears as a limit for its applicability. Proctor's 
analysis therefore does not answer what happens near this 
critical number, any more than the model in S13 does, but 
it is gratifying that it appears to have a consistent place in 
both, as well as in the numerical results. 

5.3. Dependence on Pr 

Fig. [4] shows the results of a set of simulations testing the 
dependence on Prandtl number. The Lewis number is fixed 
at Le = 0.01. The density ratio is varied over the range 
1.15 < R p < 5, Prandtl number from 0.01 to 1. Within the 
limited numerical sensitivity due to the stochastic nature of 
the flow, little systematic dependence on Pr is detectable. 
As argued above, a low dependence on Pr was to be ex- 
pected in the limit Pr 1 0. It appears that this extends 
also to a Prandtl number of order unity. 

5.4. Dependence on aspect ratio 

The aspect ratio of 2 : 1 used in the results reported above 
was tested against 5 : 1 and 10 : 1 at the same spatial res- 
olution. For the reference simulation (Pr = 0.1, Le = 0.01, 
Ra* = 10 5 , R p = 1.15) we find Nu s = 90 and Nu T = 8.5. 
By comparison the simulation with 5 : 1 results in Nus = 90 
and Nut = 9.0. The most extended box with an aspect ra- 
tio of 10 : 1 and a spatial resolution of 1500 x 300 has 
Nusselt numbers of Nus = 80 and Nut = 8.75. The aspect 
ratio thus has no significant influence on the dependences 
of the fluxes on input parameters in our simulations, within 
the fluctuations due to the stochastic nature of the flow. 
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Fig. 4. Dependence on Prandtl number and density ratio. 
Solute Nusselt number (top) and thermal Nusselt number 
(bottom) for Pr = 1, KT 1 and 10~ 2 , with Ra» = 1.6 10 5 
and Le fixed at 0.01. 



5.5. Dependence on initial stratification 

As initial state we used either a constant linear gradient 
of temperature and solute between top and bottom val- 
ues ('linear'), or a profile approximating the expected step- 
like combination of a stagnant and an overturning zone. 
The linear case shows how the oscillatory phase due to the 
Kato instability develops into overturning flow, see Figs. [5j 
[6] The duration of the initial formation process is mainly 
determined by R p and Ra*. It is of the order of 5-10 os- 
cillation periods of the initial stratification. The end state 
in both cases is statistically the same. The 'step' as initial 
condition also covers cases that are stable in linear theory 
because of the subcriticality of the system, and would not 
develop from a linear initial profile. With this initial state 
computing time can be saved for small values in Ra* and 
high R p . It was used in most of the results reported above. 

5.6. Comparison with compressible results 

The compressible simulations are based on a 5th order 
weighted ENO scheme. Compressible fluids lead to restric- 
tions in time stepping (due to the need to resolve sound 
waves). The compressible simulations take up to 100 times 
longer for the same resolution compared to the incompress- 
ible solver. In both regimes a resolution of 300 x 300 points 
has been set. Therefore only a few tests have been done 
with the compressible code. The degree of compressibility 
is governed by the e = d/H of layer thickness to pressure 
scale height. In the limit e — >• there is a direct translation 
between the compressible and the Boussinesq case (Sect. 
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Fig. 5. Example of the development of an overturning flow 
from Kato oscillations for Pr = 1, Le = 10 _1 , R p = 1.15 
and Ra* = 1.6 x 10 5 starting from a linear stratification. 
Time from left to right and top to bottom. Wave brak- 
ing occurs after 5 oscillations (Fig. 5.2). The layer is fully 
evolved after 10 oscillations (Fig. 5.4). See also the movie 
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Fig. 6. Top panel: evolution of the Nussclt number for the 
linear initial stratification. Convective cells get established 
from Kato oscillations (cf. Fig. [5]) after about 10 turnover 
times. Starting the simulation from a step (bottom) saves 
computing time. At the end of the runs the Nusselt num- 
bers of both simulations are the same within the statistical 
variations. 



6.1). The results of a numerical comparison with e = 0.1 
is shown in Table 1. The resulting Nusselt numbers do not 
differ significantly. 

Simulations done with e = 1.0 behave quite similar to 
those done with e = 0.1. The mixing processes do not sig- 
nificantly differ, at least for Pr < 0.1 and as long as the 
Rayleigh number is high enough, Ra > 5.0 x 10 5 . Differences 
of order 50% in the Nusselt numbers were found, but con- 
sidering the present level of accuracy (a factor of 2, say), we 
conclude that the Boussinesq approximation qualitatively 
gives the right results, even for layer thicknesses approach- 
ing a scale height. At Pr of order unity or larger, the effect 
of viscous damping might become more significant in more 
strongly stratified cases, however. 



Pr 


Le 


Rp 


Ra, 


Nuf 


Nu? 


Nu| 


Nut 


lO" 1 


10 -2 


2.0 


1.6 x 10 b 


60 


8 


55 


7 


lO" 1 


10" 2 


1.2 


1.6 x 10 5 


110 


12 


110 


11 


10" 1 


10" 2 


2.0 


1.6 x 10 6 


150 


12 


130 


10 


KT 1 


10~ 2 


1.2 


1.6 x 10 6 


200 


16 


200 


14 


1.0 


io- 1 


2.0 


1.6 x 10 5 


3.5 


2 


11 


1.5 


1.0 


io- 1 


1.2 


1.6 x 10 5 


45 


15 


17 


5 


1.0 


IO' 1 


2.0 


1.6 x 10 6 


4 


3.5 


26 


10 


1.0 


IO" 1 


1.2 


1.6 x 10 6 


33 


11 


26 


10 



Table 2. Comparison of compressible and incompressible 
simulations. Thermal (t) and solute (s) Nusselt numbers 
from the Boussinesq ( B ) and compressible ( s ) results. Layer 
thickness d is 0.1 pressure scale height. 



5. 7. Multi-layer simulations 

In all of the above we have assumed that the interfaces 
between the double diffusive layers can be approximated as 
solid boundaries. To test the reliability of this assumption, 
a few cases were run where the initial state consisted of two 
steps instead of a single one. In some, though not all of these 
runs, the division into two layers remained till the end. An 
example is shown in Fig. [7] for a case with Pr = 1, Le = 
10~ 2 , R p = 1.15, Ra, = 6x 10 5 and a resolution of 500 x 
500. Note the approximate (anti-)symmetry of the plumes 
near the interface in the middle, a phenomenon known from 
laboratory experiments (e.g. Fig. 1 in Turner 1985). It is 
caused by the continuity of the horizontal velocity across 
the interface enforced by viscosity. This symmetry is less 
marked in simulations at lower Prandtl numbers, when the 
thickness of the viscous boundary layer becomes smaller 
than the thickness of the stagnant zone. 

Fig. [8] shows horizontal and temporal average profiles of 
temperature and solute with height. The transition in the 
middle is broad compared with the boundary layers at top 
and bottom. Inspection of the time dependent flow shows 
that this is due to two separate effects. One is the displace- 
ments of the interface by surface waves, which smoothen the 
average gradient without changing the actual fluxes across 
the interface. In addition there is possibly some real mix- 
ing associated with breaking of the surface waves, but it 
remains localized around the interface between overturning 
and stagnant zone. 
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Fig. 7. Snapshot of a simulation including the 
free interface between two layers. Left: solute, 
Right: temperature. Pr=1.0, Le=0.01 , R p =1.15, 
Ra» = 610 5 . See movie at http://www.mpa- 



g arching . mpg.de/~henk/double Jayer . avi 



Pr=1.0 and Le=10 

Salinity 

0.6 O.t 




1.15 



0.4 0.6 
Temperature 

Fig. 8. Average temperature and solute profiles with height 
z in the double-layer simulation of Fig. [7] 

6. Application to stars 

6.1. Boussinesq limit: thin layers 

If the double-diffusive layering is thin, convective overturn- 
ing times are long compared with the sound crossing time 
of the layer thickness, and a Boussinesq (or other 'low Mach 
number' approximation) can be used for a compressible 
fluid. By the presence of a stratification in density and 
temperature, however, a heat flux is present even in a con- 
vectively neutral stratification. This requires some care in 
translating Boussinesq results to the astrophysical case. 

6.2. Heat flux 

Let F T , c be the radiative and hydrodynamic (semiconvec- 
tive) contributions to the heat flux in the star. The radiative 
heat flux is proportional to the logarithmic temperature 



gradient V = dlnT/dlnP: 

F r = fcV = £;V a + fc(V - V B ) = F IS 



(26) 



where A: is a constant depending on the local thermody- 
namic state, V a the adiabatic gradient, and F ia , F rs are 
the contributions to the radiative heat flux of the adiabatic 
and the superadiabatic parts of the mean temperature gra- 
dient, respectively. In the Boussinesq model, the contribu- 
tion F ra is absent: the convective and radiative heat fluxes 
F c , F r are governed by the same temperature gradient. 
Related to this, the Boussinesq model has one parame- 
ter less: the pressure scale height H. Because of this dif- 
ference, the heat flux in the stellar (compressible) model 
cannot be compared directly with the heat flux in an in- 
compressible model. Instead, the ratio of convective to ra- 
diative heat flux in the Boussinesq model is to be identified 
with the ratio / = F c /F Ta of the convective flux to the su- 
peradiabatic component of the radiative flux in the star, 
in the limit H -> oo (cf. Massaguer & Zahn 1980). The 
Boussinesq model thus is the limit e | (taken at fixed 
Ra*). If the semiconvective layer thickness e is 'sufficiently 
small', the compressible case can be compared directly with 
the Boussinesq model (as veri fied by the numerical tests 
with e = 1 and e = 0.1 in Sect. 5.6 above). 



This makes the semiconvection problem more amenable 
to numerical simulation. In contrast with a convective stel- 
lar envelope for example, with its many scale heights to be 
covered, the layered nature of double diffusive convection 
puts it in a parameter range that is much closer to condi- 
tions accessible with realistic ab initio calculations. 



With ( 26 ) for the radiative flux, the total heat flux can 



be written as 



F = F T + F conv = fcV a + Nu fc(V 



(27) 



where Nu is the Nusselt number. By taking out the radia- 
tive flux associated with the adiabatic gradient, the Nusselt 
number as defined in (27 1 can now be identified with the 
Boussinesq equivalent (see also Massaguer & Zahn 1980). 
The relevant Nusselt Number is thus not simply the ratio 
of total to radiative heat fluj|3 

6.3. Rayleigh number and density ratio 

The modified Rayleigh number for a layer of thickness d is, 
in astrophysical notation: 



Ra* = Pr Ra 



gd 4 

n 2 H 



(V- V, 



(28) 



where k is the thermal diffusivity (cm 2 /s) and g the acceler- 
ation of gravity. Ra* is typically a very large number unless 
the layer thickness d is small (for estimates see the 15 M 
example below). Apart from Le, the problem also depends 
on the relative strength of the stabilizing solute gradient 
relative to the destabilizing thermal gradient. This can be 
measured in terms of the thermal and solute buoyancy fre- 
quencies TVt, N$: 



fJ 



JVt = ij(V.-V), 



H 



(29) 
(30) 



2 The implications of this distinction are not apparent in some 
of the published work on semiconvection. 
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with 



dln/x/dlnP, 



(31) 



where fi is the mean atomic weight per particle, and the 
gradients are understood as mean values over a double dif- 
fusive step. The density ratio R p : 



R„ 



(32) 



is then the equivalent of the density ratio in the Boussinesq 
formulation, and is the dimensionless measure we will use 
for the relative strength of the stabilizing Helium gradient. 
Semiconvection, i.e. a stratification 'between Schwarzschild 
and Ledoux' then corresponds to R p > 1. V— V a is typically 
a small number in a stellar interior (even in the presence 
of the double diffusive steps) , and R p a large number (see 
below) . 

6.4. Extrapolation 

The connection between the numerical results and the as- 
trophysical conditions requires extrapolation that cannot at 
present be covered by the simulations. It is covered by the 
model in S13, as validated by the comparison with our nu- 
merical results. In contrast with the laboratory setup where 
the temperature difference is the imposed quantity, the heat 
flux is the fixed quantity under the conditions in a stellar 
semiconvective zone, since the structure of the star depends 
only little on the behavior of the semiconvective zone. The 
asymptotic conditions for large Ra* then lead (S13 Sect. 6) 
to simple expressions for the superadiabaticity V — V a and 
the effective Helium diffusivity K se g: 



V - V a « Le 1/2 V, 



K scff — ( K S K 



,1/2 (Vr- Va) 



(33) 
(34) 



(S13 eqs. 57, 59). The analysis in S13 does not take into ac- 
count a modification due to the effect of r adia tion pressure 
discussed in S92. With this modification, (34 1 becomes: 



Kscff = (k s k) 1/2 (- - 3) 



(V r 



(35) 



where j3 = P g / (P g + P r ) is the ratio of gas to total pressure. 
The density ratio approaches the value 



R„ 



Le- 1 / 2 , 



(36) 



while the relative thickness of the stagnant zones in this 
limit is of order (setting (3 — 1) 



S 1/Nu s = Le 1/2 V M /(V r 
The range of validity of the limit is 
l < d < H, 

where 



(37) 



(38) 



l = (k 2 H/ 9 )V 4 (39) 

is the length scale on which the thermal diffusion time scale 
equals the free fall time over a pressure scale height. 

The asymptotic value R p = Le -1 / 2 (eq. 36) is also the 
maximum density ratio for which the theory predicts exis- 
tence of the layered state. Semiconvection in a stellar model 



with heat flux fixed is thus predicted to settle close to this 
maximum, if the layer thickness is not small (d 3> do or 
Ra* ^> 1). How close it actually settles can also be pre- 
dicted, but requires analysis to next order in the small 
quantity 1/Ra*. 

In a stellar evolution model, expressions p3|35l ) give 
nearly identical results as the model of Zaussinger (2011), 
which is based on S92. They are easier to implement since 
they have been restricted to conditions d ^> lo that are the 
most interesting for stellar evolution anyway. 



6.5. A 15 M® star 

We are now in the position to estimate the range of pa- 
rameter values for semiconvection in a star. Consider the 
important case of massive stars around main sequence 
turnoff. We use a model kindly provided by A. Weiss (model 
'fzml5_151'). Characteristic values for the physical quanti- 
ties in the semiconvective zone of this model are q m 10 6 



cm/s , Kg 1 cm 2 /s, k ps 3 • 10 8 cm 2 /s, H ps 2 
V a = 0.4, V r - V a « 0.02, V, 



9 
10 10 



1. 



With ( 28 1 and ( 33 ) , the Rayleigh number can be found 



for this semiconvective zone, as a function of the layer thick- 
ness d. This yields Ra* - 10 12 for d/H = 0.1, or Ra, - 10 8 
for d/H = 0.01, for example. For such 'macroscopic' layer 
thicknesses, conditions are thus in the asymptotic regime 
for which the expressions in 6.4 are approximately valid. 
The effective He-diffusivity from (35) is k sc s ~ 10 3 



cm /s, three orders of magnitude above the microscopic 
value. The mixing time scale over a pressure scale height is 
thus about 10 10 yr. The value of is 2 10 5 cm, or ab out 
10~ 5 pressure scale heights. The limiting expressions (33 
37 1 are thus applicable over a large range in the value of the 
(uncertain) layer thickness. The lower end of this range is 
not likely to be very relevant, since at d ~ Iq the expected 
time scale for layer merging (see below) is very short, of the 
order of Le -1 / 2 times the free fall time scale over a pressure 
scale height (from Eqs. 39 42). 

The time scale for the initial layer formation from a 
Kato oscillation is similarly short, a few times the instabil- 
ity growth time, which is of the order of the convective over- 
turning time in the absence of the stabilizing /^-gradient. 
For the 15 M star this works out to about one day. The 
asymptotic expressions ( 33 1 and ( 35 1 should therefore be 



adequate for use in stellar evolution codes 



6.6. Layer thickness 

While the layer thickness has little influence on the effective 
mixing rate, it is useful to check if it could become signif- 
icant compared with the pressure scale height, in which 
case the limit of thin layering assumed would not be valid. 
This requires a model for the layer thickness, for which no 
good theory is available. Observations in laboratory exper- 
iments and in geophysical cases show that layer thickness is 
not constant, but grows in time by processes of merging of 
neighboring layers (McDougall 1981, Ross & Lavery 2009, 
see also the numerical experiment in Young & Rosner 2000). 
Layer thickness can therefore not be treated independently 
of the history of the system. 

An estimate can be made, however, from the effective 
solute diffusivity. Changes in layer thickness involve the ex- 
change of solute between neighboring layers. Over a dis- 
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tance D, the time scale r for exchange is given by the effec- 
tive diffusivity K ae s: t m _D 2 /re se ff. Setting D equal to the 
layer thickness d itself then gives an estimate of the rate of 
change of thickness by merging: 



or 



dlnd/dt = 1/r, 



d , Kgeff 
— d = — . 

d* d 



If k s e ff is constant in time: 

d = (2 KscS t) 1 / 2 , 



(40) 



(41) 



(42) 



where t is time since the semiconvective condition started. 
In the 15 M star the duration of the semiconvective state is 
of the order of 10 7 yr, so the layer thickness to be expected 
at the end of the semiconvective phase is of the order 10 9 
cm, or 0.05 H. Given the uncertainties involved, it cannot 
be ruled out that semiconvective layer thickness d can ac- 
tually become 'macroscopic', of order H, in the course of 
the star's evolution. 



7. Discussion 

As in the case of ordinary convection, there are far fewer 
intrinsic parameters in semiconvection than the quantities 
defining the physical state in a stellar interior. This allows a 
significant volume of astrophysically useful parameter space 
to be covered by a grid of numerical simulations. A physical 
model as applied here is still needed, however, to interpret 
the results and to extrapolate them meaningfully to the 
astrophysically relevant range. 

The low value of the viscosity and solute diffusivity com- 
pared with the thermal diffusivity constitute a limiting case 
that actually simplifies the double diffusive problem greatly. 
Among other things, the results become nearly independent 
of the Prandtl number in this limit [as suggested already 
by Proctor's (1981) analysis]. 

The simulations were all done in 2-D, so tests in 3-D 
will be needed for verification. It is unlikely, however, that 
the results will turn out very different, at least within an 
astrophysical factor of 2. The reason is that in the low-Pr, 
low-Le limit the flow in the layers is almost equivalent to 
ordinary convection between plates. Due to the low solute 
diffusivity, the amount of solute in the bulk of the layer is 
small. It can then be treated as a perturbation, as assumed 
in S92. Known laboratory results for convection at very 
high Rayleigh numbers can then be used to extrapolate the 
numerical results. As shown in S13, this makes predictions 
similar to the simple 2-D model used in S92. In particular 
the effective mixing rate is very low, as also found (in a dif- 
ferent region in parameter space) in geophysical examples 
like lake Kivu and the (ant)arctic oceans (cf. Sect. [2] and 
S13). 

Stochastic fluctuations in the flow produce scatter in 
the fluxes of heat and solute, which affect the measured 
averages (cf. Table 1). The accuracy of these averages could 
be improved with longer runs. 

Most of the results presented are based on simulation of 
a single double diffusive layer. As we found in |5.7| this does 
not fully reproduce broadening of the interfaces between 
layers by surface waves (see movie at Fig. \fh . 

The physics determining the thickness of the individual 
double diffusive layers remains uncertain. In the equivalent 



geophysical examples semiconvective zones always consist 
of a large number of very long-lived layers. Observations 
in the east-african rift lakes (e.g. Schmid et al. 2010) show 
that layers first forming at the boundary of an expanding 
double-diffusive zone are always thin, subsequently growing 
slowly by a process of merging. This is understood in terms 
of an energy argument (cf. Sect. [2]), and is likely to happen 
in a growing stellar se mico nvective zone as well. 

The estimate (eq. [42] ) suggests that layer thickness 
might approach a significant fraction of a pressure scale 
height. In this case, only a few layers would be left at 
the end of the semiconvective phase, and the location of 
their boundaries may well vary somewhat randomly be- 
tween stars. As surmised in Zaussinger (2011), this might 
introduce a random element in the late stages of the evolu- 
tion of massive stars. 

Semiconvection as studied here addresses only one of 
two astrophysical meanings of the term. In addition to the 
effect of a stabilizing solute gradient, there is the effect of 
a composition-dependent opacity: the increase of the ra- 
diative gradient when mixing is assumed, with consequent 
onset of convection in a stratification which was radiative 
before mixing. Historically, this has actually been the main 
concern in stellar evolution computations. The physics of 
this kind of semiconvection is completely different from the 
double-diffusive kind (cf. Kippenhahn et al. 2013); it has 
not received the same theoretical attention. 

8. Conclusions 

The results show how the physics of double diffusive convec- 
tion known from geophysical and laboratory studies can be 
applied to the astrophysical case of low Prandtl and Lewis 
numbers. As expected, the process takes place in the form 
of the characteristic double diffusive layering known and 
theoretically understood since the 1980 's. 

In the asymptotic regime occupied by astrophysical con- 
ditions the number of independent parameters determining 
the physics is effectively reduced to three, so that a mean- 
ingful range of parameter values can be covered with numer- 
ical simulations. We have compared a grid of simulations 
with the predictions of the model for a semiconvective layer 
in S13. In the parameter region where the estimate over- 
laps with the region covered by the simulations we find 
good agreement, even without significant tuning of the one 
fitting parameter used (the 'erosion factor' q, Sect. 3.5 1. 

The simulations also provide an answer to an interesting 
question raised by Proctor's (1981) mathematical analysis 
and by the model in S13. Both predicted that a critical 



number R n 



Le 



-1/2 



plays a role in the steady layered 



state. The results presented in Sect. 5.2 clarify the behavior 
of the layered state around this maximum density ratio. In 



a system set up in a layered state just above R„ 



the 



Nusselt numbers show large fluctuations in amplitude; the 
duration of the fluctuations becomes increasingly long as 
R p approaches R pmax - This explains why theories based on 
the assumption of stationarity fail near this limit. 

Boussinesq and fully compressible simulations give 
equivalent results in the limit of small layer thickness. The 
Boussinesq calculations even reproduce the approximate 
Nusselt numbers for layers as thick as a scale height. The 
numerical results confirm the theoretical expectation that 
at low Prandtl number the dependence on viscosity is weak, 
and extends to Pr of order unity. 
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The effective mixing rate in the semiconvective zone of a 
15 Mq main sequence star is predicted to be low, though the 
semiconvective phase may possibly leave significant jumps 
in the profile of Helium concentration. 
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